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Abstract The radioscience experiment is one of the on board experiment of the Mer- 
cury ESA mission BepiColombo that will be launched in 2014. The goals of the experi- 
ment are to determine the gravity field of Mercury and its rotation state, to determine 
the orbit of Mercury, to constrain the possible theories of gravitation (for example 
by determining the post-Newtonian (PN) parameters), to provide the spacecraft posi- 
tion for geodesy experiments and to contribute to planetary ephemerides improvement. 
This is possible thanks to a new technology which allows to reach great accuracies in 
the observables range and range rate; it is well known that a similar level of accuracy 
requires studying a suitable model taking into account numerous relativistic effects. 
In this paper we deal with the modelling of the space-time coordinate transformations 
needed for the light-time computations and the numerical methods adopted to avoid 
rounding-off errors in such computations. 

Keywords Mercury • Interplanetary tracking • Light-time • Relativistic effects • 
Numerical methods 



1 Introduction 

BepiColombo is an European Space Agency mission to be launched in 2014, with the 
goal of an in-depth exploration of the planet Mercury; it has been identified as one 
of the most challenging long-term planetary projects. Only two NASA missions had 
Mercury as target in the past, the Mariner 10, which flew by three times in 1974-5 and 
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Messenger, which carried out its flybys on January and October 2008, September 2009 
before it starts its year-long orbiter phase in March 2011. 

The BepiColombo mission is composed by two spacecraft to be put in orbit around 
Mercury. The radioscience experiment is one of the on board experiments, which would 
coordinate a gravimetry, a rotation and a relativity experiment, using a very accurate 
range and range rate tracking. These measurements will be performed by a full 5- 
way link (0]) to the Mercury orbiter; by exploiting the frequency dependence of the 
refraction index, the differences between the Doppler measurements (done in Ka and X 
band) and the delay give information on the plasma content along the radiowave path. 
In this way most of the measurements errors introduced can be removed, improving of 
about two orders of magnitude with respect to the past technologies. The accuracies 
that can be achieved are 10 cm in range and 3 x 10~ 4 cm/s in range rate. 

How we compute these observables? For example, a first approximation of the range 
could be given by the formula 

r = |r| = |(x sat + x M ) - (x EM + x E + x ant )| , (1) 

which models a very simple geometrical situation (see Figure [1]). The vector x sa t is the 
mercurycentric position of the orbiter, the vector is the position of the center of 
mass of Mercury (M) in a reference system with origin at the Solar System Barycenter 
(SSB), the vector x^m is the position of the Earth- Moon center of mass in the same 
reference system, xg is the vector from the Earth-Moon Barycenter (EMB) to the 
center of mass of the Earth (E) , the vector x an t is the position of the reference point 
of the ground antenna with respect to the center of mass of the Earth. 




SSB 



Fig. 1 Geometric sketch of the vectors involved in the computation of the range. SSB is the 
Solar System Barycenter, M is the center of Mercury, EMB is the Earth-Moon Barycenter, E 
is the center of the Earth. 

Using ([l| means to model the space as flat arena (r is an Euclidean distance) and 
the time as absolute parameter. This is obviously not possible because it is clear that, 
beyond some threshold of accuracy, these quantities have to be formulated within the 
framework of Einstein's theory of gravity (general relativity theory, GRT) . Moreover we 
have to take into account the different times at which the events have to be computed: 
the transmission of the signal at the transmit time (it), the signal at the Mercury 
orbiter at the time of bounce (tf,) and the reception of the signal at the receive time 
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Formula {TJ could be a good starting point to construct a correct relativistic for- 
mulation; with the word "correct" we do not mean all the possible relativistic effects, 
but the effects that can be measured by the experiment. This paper deals with the 
corrections to apply to this formula to obtain a consistent relativistic model for the 
computations of the observables and the practical implementation of such computa- 
tions. 

In Section [2] we discuss the relativistic four-dimensional reference systems used 
and the transformations adopted to make the sums in fT} consistent; according to 
[l^ . with "reference system" we mean a purely mathematical construction, while a 
"reference frame" is a some physical realization of a reference system. The relativistic 
contribution to the time delay due to the Sun's gravitational field, the Shapiro effect, 
is described in Section [3] Section [4] deals with the theoretical procedure to compute 
the light-time (range) and the Doppler shift (range rate). In Section [5] we discuss 
the practical implementation of the algorithms showing how we eliminate rounding-off 
problems. 

The equations of motion for the planets Mercury and Earth, including all the 
relativistic effects (and potential violations of GRT) required to the accuracy of the 
BepiColombo radioscience experiment have already been discussed in [Io| . thus this 
paper concentrates on the computation of the observables. 

2 Space-time reference frames and transformations 

The five vectors involved in formula fl} have to be computed at their own time, the 
epoch of different events: e.g., x an t, x^m an d xg are computed at both the antenna 
transmit time tt and the receive time t r of the signal, xjyj and x sa t are computed at 
the bounce time t^ (when the signal has arrived to the orbiter and is sent back, with 
correction for the delay of the transponder). To be able to perform the vector sums 
and differences, these vectors have to be converted to a common space-time reference 
system, the only possible choice being some realization of the BCRS (Barycentric Ce- 
lestial Referece System). We adopt for now a realization of the BCRS that we call SSB 
(Solar System Barycentric) reference frame and in which the time is TDB (Barycentric 
Dynamic Time); other possible choices, such as a TCB (Barycentric Celestial Time), 
only can differ by linear scaling. The TDB choice of the SSB timescale entails also the 
appropriate linear scaling of space-coordinates and planetary masses as described for 
instance in 0] or 0]. 

The vectors xy, xg, and xjm are already in SSB as provided by numerical integra- 
tion and external ephemerides; thus the vectors x an t and x sa t have to be converted to 
SSB from the geocentric and mercurycentric systems, respectively. Of course the con- 
version of reference system implies also the conversion of the time coordinate. There are 
three different time coordinates to be considered. The currently published planetary 
ephemerides are provided in TDB. The observations are based on averages of clock 
and frequency measurements on the Earth surface: this defines another time coordi- 
nate called TT (Terrestrial Time) . Thus for each observation the times of transmission 
tt and receiving t r need to be converted from TT to TDB to find the corresponding 
positions of the planets, e.g., the Earth and the Moon, by combining information from 
the precomputed ephemerides and the output of the numerical integration for Mercury 
and the Earth-Moon barycenter. This time conversion step is necessary for the accurate 
processing of each set of interplanetary tracking data; the main term in the difference 
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TT-TDB is periodic, with period 1 year and amplitude ~ 1.6 x 10~ 3 s, while there is 
essentially no linear trend, as a result of a suitable definition of the TDB. 

The equation of motion of a mercurycentric orbiter can be approximated, to the 
required level of accuracy, by a Newtonian equation provided the independent variable 
is the proper time of Mercury. Thus, for the BepiColombo radioscience experiment, 
it is necessary to define a new time coordinate TDM (Mercury Dynamic Time), as 
described in [10|], containing terms of 1-PN order depending mostly upon the distance 
from the Sun and velocity of Mercury. 

From now on we shall call the quantities related to the SSB frame "TDB-compatible" , 
the quantities related to the geocentric frame "TT-compatible" , and the quantities re- 
lated to the mercurycentric frame "TDM-compatible", in accordance with the paper 
7], and label them TB, TT and TM, respectively. 

The differential equation giving the local time T as a function of the SSB time t , 
which we are currently assuming to be TDB, is the following: 



— - \ — 

~dt ~ c2 



v 2 

U+--L 



(2) 



where U is the gravitational potential (the list of contributing bodies depends upon 
the accuracy required: in our implementation we use Sun, Mercury to Neptune, Moon) 
at the planet center and v is the SSB velocity of the same planet. The constant term 
L is used to perform the conventional rescaling motivated by removal of secular terms, 
e.g., for the Earth we use Lq. 

The space-time transformations we have to perform involve essentially the position 
of the antenna and the position of the orbiter. The geocentric coordinates of the an- 
tenna should be transformed into TDB-compatible coordinates; the transformation is 
expressed by the formula 

/ TT \ 

TB TT U , \ 1 / V E • X ant \ TB 



where U is the gravitational potential at the geocenter (excluding the Earth mass), 
Lq = 1.48082686741 x 10~ 8 is a scaling factor given as definition, supposed to be 
a good approximation for removing secular terms from the transformation and v™ 
is the barycentric velocity of the Earth. The next formula contains the effect on the 
velocities of the time coordinate change, which should be consistently used together 
with the coordinate change: 
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Note that the previous formula contains the factor dT/dt (expressed by eq. @) that 
deals with time transformation: T is the local time for Earth, that is TT, and t is the 
corresponding TDB time. 

The mercurycentric coordinates of the orbiter should be transformed into TDB- 
compatible coordinates through the formula 

TB _ TM /, U \ 1 / Vm B ' X-Tat 1 \ TB 
X S at — X sat I 1 - ^ - LCM I - g I ^2 ) M ' 
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where U is the gravitational potential at the center of mass of Mercury (excluding the 
Mercury mass) and Lqm could be used to remove the secular term in the time trans- 
formation (thus defining a TM scale, implying a rescaling of the mass of Mercury). We 
believe this is not necessary: the secular drift of TDM with respect to other time scales 
is significant, see Figure 5 in paper [Io| . but a simple iterative scheme is very efficient 
in providing the inverse time transformation. Thus we set -LcAf = 0, assuming the 
reference frame is TDM-compatible. As for the antenna we have a formula expressing 
the velocity transformation that contains the derivative of time T for Mercury, that is 
TDM, with respect to time t, that is TDB: 



TB _ 
v sat — 



TM 
v sat 
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1 2 — L CM 
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In all the formulas for these coordinate changes we have neglected the terms of 
the SSB acceleration of the planet center (01), because they contain beside 1/c 2 the 
additional small parameter (distance from planet center) /(planet distance to the Sun), 
which is of the order of 10~ 4 even for a Mercury orbiter. 

To assess the relevance of the relativistic corrections of this section to the accuracy 
of the BepiColombo radioscience experiment, we have computed the observables range 
and range rate with and without these corrections. As shown in Figure(2] the differences 
are significant, at a signal-to-noise ratio S/N ~ 1 for range, much more for range rate, 
with an especially strong signature from the orbital velocity of the mercurycentric orbit 
(with S/N > 50). 



3 Shapiro effect 

The correct modelling of space-time transformations is not sufficient to have a precise 
computation of the signal delay: we have to take into account the general relativistic 
contribution to the time delay due to the space-time curvature under the effect of the 
Sun's gravitational field, the Shapiro effect ([I2])- The Shapiro time delay At at the 
1-PN level, according to [l5j] and [ll]], is 



At = ( 1+ 7^° 



rt — |rt| and r r = |r r | are the heliocentric distances of the transmitter and the receiver 
at the corresponding time instants of photon transmission and reception, fio is the grav- 
itational mass of the Sun (/j,q = Gmo) and r = |iy — rt|. The planetary terms, similar 
to the solar one, can also be included but they are smaller than the accuracy needed 
for our measurements. Parameter 7 is the only post-Newtonian parameter used for the 
light-time effect and, in fact, it could be best constraint during superior conjunction 
([3|). The total amount of the Shapiro effect in range is shown in Figure [3] 

The question arises whether the very high signal to noise in the range requires other 
terms in the solar gravity influence, due to either (i) motion of the source, or (ii) higher- 
order corrections when the radio waves are passing near the Sun, at just a few solar 
radii (and thus the denominator in the log- function of the Shapiro formula is small). 
The corrections (i) are of the post-Newtonian order 1.5, that is containing a factor 
I/c 3 , but it has been shown in [l^] that they are too small to affect our accuracy. 
The corrections (ii) are of order 2, that is 1/c 4 , but they can be actually larger for 
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Fig. 2 The difference in the observables range and range rate for one pass of Mercury above 
the horizon for a ground station, by using an hybrid model in which the position and velocity of 
the orbiter have not transformed to TDB-compatible quantities and a correct model in which 
all quantities are TDB-compatible. Interruptions of the signal are due to spacecraft passage 
behind Mercury as seen for the Earth station. Top: for an hybrid model with the satellite 
position and velocity not transformed to TDB-compatible. Bottom: for an hybrid model with 
the position and velocity of the antenna not transformed to TDB-compatible. 



an experiment involving Mercury. The relevant correction is most easily obtained by 
adding 1/c 4 terms in the Shapiro formula, due to the bending of the light path: 

rt + r? . + r+ (i±7)w \ 

This formulation has been proposed in [llj] and it has been justified in the small 
impact parameter regime by much more theoretically rooted derivations in [13 ] 
and [J. Figure [4] shows that the order 2 correction is relevant for our experiment, 
especially when there is a superior conjunction with a small impact parameter of the 
radio wave path passing near the Sun. Note that the 1/c 4 correction (~ 10 cm) in the 



S(i) 



(1 + 7) MO 



In 



7 



tg 6 Change in the observable 




100 200 300 400 500 600 700 800 

time, days from arc beginning 



Fig. 3 Total amount of the Shapiro effect in range over 2-year simulation. The sharp peaks 
correspond to superior conjunctions, when Mercury is "behind the Sun" as seen from Earth, 
with values as large as 24 km for radiowaves passing at 3 solar radii from the center of the 
Sun. Interruptions of the signal are due to spacecraft visibility from the Earth station (in this 
simulation we assume just one station). 



Shapiro formula effectively corresponds to ~ 3 x 10~ J correction in the value of the 
post-Newtonian parameter 7. 

The Shapiro correction for the computation of the range rate is 



S = 



2(1+ 7 )W) 



-r {ft + f r ) + f (r t + r r + 



where the time derivative is with respect to a TDB time. This formula is almost never 
found in the literature and has not been much used in the processing of the past 
radioscience experiments, such as in Q], because the observable range rate is typically 
computed as difference of ranges divided by time; however, for reasons explained in 
Section [5] this formula is now necessary. 



4 Light-time iterations 

Since radar measurements are usually referred to the receive time t r the observables are 
seen as functions of this time, and the computation sequence works backward in time: 
starting from t r , the bounce time is computed iteratively, and, using this information 
the transmit time tt is computed. 

rp 73 rp rj 

The vectors xj^ and xg^ are obtained integrating the post-Newtonian equations 
of motions. The vectors x^f are obtained by integrating the orbit in the mercurycen- 
tric TDM-compatible frame. The vector x^ is obtained from a standard IERS model 
of Earth rotation, given accurate station coordinates, and Xg T from lunar ephemerides 

(0)- 

In the following subsections we shall describe the procedure to compute the range 
(Section [4T} and the range rate (Section I4.2[) . 
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Fig. 4 Differences in range (top) and range rate (bottom) by using an order 1 and an order 
2 post-Newtonian formulation (7 = 1); the correction is relevant for BepiColombo, at least 
when a superior conjunction results in a small impact parameter b. E.g., in this figure we 
have plotted data assumed to be available down to b ~ 3Ho- For larger values of b the effect 
decreases as 1/b 2 . 



4.1 Range 

Once the five vectors are available at the appropriate times and in a consistent SSB 
system, there are two different light-times, the up-leg At U p — t^—tt for the signal from 
the antenna to the orbiter, and the down-leg At down — t r — t b for the return signal. 
They are defined implicitly by the distances down-leg and up-leg 

Tdo{U) = Xsat(*6(*r)) + *M(tb{tr)) - x E m(M - x B (t r ) - x ant (i r ) , 

rdo{tr) = \r do {t r )\ , c(t r - t b ) = r do (t r ) + S do (7) , (3) 



r up (t r ) = x sat (i b (i r )) + x M (t b (i r )) - x EM (t t (i r )) - x E (t t (t r )) - Xant (*t(*r)) , 
r up (t r ) = \r up {t r )\ , c(t b - t t ) = r up (t r ) + Supil) , (4) 
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Fig. 5 The difference in the observables range and range rate using a light-time in TT and a 
light time in TDB: the difference in range is very high, more than 13 meters in one day, while 
the difference in range rate is less than the accuracy of the experiment. 



respectively, with somewhat different Shapiro effects. Then tr — tf, and t^ — tt axe the 
two portions of the light-time, in the time attached to the SSB, that is TDB; this 
provides the computation of tt. Then these times are to be converted back in the time 
system applicable at the receiving station, where the time measurement is performed, 
which is TT (or some other form of local time, such as the standard UTC) . t r is already 
available in the local time scale, from the original measurement, while tt needs to be 
converted back from TDB to TT. The difference between these two TT times is Attot, 
from which we can conventionally define r(t r ) = cAttot/2- Note that the difference 
Attot in TT is significantly different from tr ~ tt in TDB, by an amount of the order 
of 10 -7 s, while the sensitivity of the BepiColombo radioscience experiment is of the 
order of 10~ 9 s, thus these conversions change the computed observable in a significant 
way, see Figure [5] 

The practical method for solving tf,(t r ) and tt(t r ) in eqs. Q and @ is as follows. 
Since the measurement is labeled with the receive time t r , the iterative procedure 
needs to start from eq. (|3]) by computing the states xem> x E an d x an t at epoch t r , 
then selecting a rough guess t9 for the bounce time (e.g., t® — t r ). Then the states x sa t 
and xm are computed at t® and a successive guess t\ is given by (J3J). This is repeated 
computing t%, and so on until convergence, that is, until t^ — is smaller than the 
required accuracy. This fixed point iteration to solve the implicit equation for tf, is 
convergent because the motion of the satellite and of Mercury, in the time t r —tf,, is 
a small fraction of the total difference vector. After accepting the last value of t^ we 
start with the states x sa t and xjj at % and with a rough guess ij for the transmit time 
(e.g., t^ = tf)). Then xem> x e an d x an t are computed at epoch t^ and t\ is given by 
eq. (d}, and the same procedure is iterated to convergence, that is to achieve a small 
enough tt — 1 . This double iterative procedure to compute range is consistent with 
what has been used for a long time in planetary radar, as described in the paper [l6| . 
We conventionally define r = (r do + S^o + 



10 



4.2 Range rate 



After the two iterations providing at convergence and tt are complete, we can proceed 
to compute the range rate. We rewrite the expression for the Euclidean range (down-leg 
and up- leg) as a scalar product: 



2 u 
r up {tr 



[X-Msih) -X Ea (t r )] ■ [x Ms (i(,)-X Ea (tr)] 
[x Ms {t b ) - X Ea (tt)] ■ [xMsOfc) - X Ea (*t)] . 



where xjj s = + x sa t and x Ea = x E m + x E + x all t . The light-time equation contains 
also the Shapiro terms, thus the range rate observable contains also additive terms S do 
and S U p, with significant effects (a few cm/s during superior conjunctions). Since the 
equations giving t b and tt are still ([3} and (O, in computing the time derivatives, we 
need to take into account that t b = tb{t r ) and tt = ii(t-r), with non-unit derivatives. 

Computing the derivative with respect to the receive time t r , and using the dot 
notation to stand for d/dt r , we obtain: 



fdo(tr) = ido 



' up \ 



i up 



x Ms(ifc 



x Ms(ifc 



1 - 



fdo(tr) + S do 



- x Ea (t r 



(5) 



fdo(tr) + S do 



X Ea (tt) 1 



rdoitr) + S do r up (t r ) + Sup 



(6) 



where 



Tdo 



x Ms(*b) - X Ea (t r ) 



XMs(*&) - X Ea (tt) 



^do(ir) " "' r U p(tr) 

However, the contribution of the time derivatives of the Shapiro effect to the d%jdt r 
and dtt/dtr corrective factors is small, of the order of 10~ 10 , which is marginally 
significant for the BepiColombo radioscience experiment. We conventionally define r = 
c(l — i{)/2 = (r do + Sdo + fup + S U p)/2. These equations are compatible with 
taking into account that they use a single iteration. 

Since the time derivatives of the Shapiro effects contain and r r , the equations 
([5]) and © are implicit, thus we can again use a fixed point iteration. It is also possible 
to use a very good approximation which solves explicitly for r do and then for r U p, 
neglecting the very small contribution of Shapiro terms: 



Tdo = r do 



x Ms( i b) 



) - X L ; ;1 (f, 



XMs(^fc) ' ido 



where the right hand side is weakly dependent upon r do only through S do , thus a 
moderately accurate approximation could be used in the computation of S do , followed 
by a single iteration. For the other leg 
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All the above computations are in SSB with TDB; however, the frequency mea- 
surements, at both tt and tr, are done on Earth, that is with a time which is TT. This 
introduces a change in the measured frequencies at both ends, and because this change 
is not the same (the Earth having moved by about 3 x 10 -4 of its orbit) there is a 
correction needed to be performed. The quantity we are measuring is essentially the 
derivative of tt with respect to t r , but this in two different time systems: for readability, 
we use T for TT, t for TDB 



dT t 
dT r 



m 

dt t 



dk 

dt r 



dt r 
dT r 



where the derivatives of the time coordinate changes are the same as the right hand 
sides of the differential equation giving T as a function of t in the first factor and the 
inverse of the same for the last factor. However, the accuracy required is such that the 
main term with the gravitational mass of the Sun hq and the position of the Sun xo is 
enough: 



dT t 
dT r 



/'n 



l*E(it)| 2 



\x E (tt)-xo{tt)\c 2 
Mo 

|x E (*r) ~ X (tr)| c 2 



2 c 2 

jx E (Mi 2 

2 c 2 



dtt_ 

dt r 
-1 



(7) 



Note that we do not need the Lq constant term discussed above because it cancels in 
the first and last terms in the right hand sides of eq. (J7J). The correction in the above 
formula is required for consistency, but in fact the correction has an order of magnitude 
of 10 -7 cm/s and is negligible for the sensitivity of the BepiColombo radioscience 
experiment (Figure [5]). 



5 Numerical problems and solutions 

The computation of the observables, as presented in the previous section, is already 
complex, but still the list of subtle technicalities is not complete. 

A problem well known in radioscience is that, for top accuracy, the range rate mea- 
surement cannot be the instantaneous value r(t r ) = (^do(^r) + S^o + r U p(t r ) + S U p)/2. 
In fact, the measurement is not instantaneous: an accurate measure of a Doppler effect 
requires to fit the difference of phase between carrier waves, the one generated at the 
station and the one returned from space, accumulated over some integration time A, 
typically between 10 and 1000 s. Thus the observable is really a difference of ranges 

• x r(t b + A/2) - r(t b - A/2) 
or, equivalently, an averaged value of range rate over the integration interval 



rt b +A/2 

TA{tr) = -r / r{s)ds. (9) 

^ >t b -A/2 



In order to understand the computational difficulty we need to take also into ac- 
count the orders of magnitude. As said in the introduction, for state of the art tracking 
systems, such as those using a multi-frequency link in the X and Ka bands, the accu- 
racy of the range measurements can be ~ 10 cm and the one of range-rate 3 x 10 -4 
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cm/s (over an integration time of 1000 s). Let us take an integration time A = 30 
s, which is adequate for measuring the gravity field of Mercury; in fact if the orbital 
period is ~ 8 000 s, the harmonics of order m — 26 have periods as short as ~ 150 s. 
The accuracy over 30 s of the range rate measurement can be, by Gaussian statistics, 
~ 3 x 10~ 4 ^1000/30 ~ 17 x 10~ 4 cm/s, and the required accuracy in the compu- 
tation of the difference r(t{, + A/2) — r(i(, — A/2) is ~ 0.05 cm. The distances can 
be as large as ~ 2 x 10 13 cm, thus the relative accuracy in the difference needs to be 
2.5 x 10~ lj . This implies that rounding off is a problem with current computers, with 



relative rounding off error of e — 2 



-o2 



2.2 x 10 



-16 



(Figure extended precision 



is supported in software, but it has many limitations. The practical consequences are 
that the computer program processing the tracking observables, at this level of preci- 
sion and over interplanetary distances, needs to be a mixture of ordinary and extended 
precision variables. Any imperfection may result in "banding", that is residuals show- 
ing a discrete set of values, implying that some information corresponding to the real 
accuracy of the measurements has been lost in the digital processing. 



Charge in the observable 



}, seconds from arc beginning 
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Fig. 6 Range and range rate differences due to a change by 10" 11 of the C'22 harmonic 
coefficient: the range rate computed as range difference divided by the integration time of 30 
s, eq. is obscured by the rounding off. 



As an alternative, the use of a quadrature formula for the integral in eq. © can 
provide a numerically more stable result, because the S/N of the range rate measure- 
ment is <C 1/e. Figure [7] shows that a very small model change, generating a range 
rate signal < 2 micron/s over one pass, can be computed smoothly by using a 7 nodes 
Gauss quadrature formula. 



6 Conclusions 



By combining the results of the previous paper ) 5 an d of this one, we have com- 
pleted the task of showing that it is possible to build a consistent relativistic model of 
the dynamics and of the observations for a Mercury orbiter tracked from the Earth, at 
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Change in the observable 



}, seconds from arc beginning 




time, seconds trom arc beginning 



Fig. 7 Range and range rate differences due to a change by 10 -11 of the C22 harmonic coef- 
ficient: the range rate computed as an integral, eq. |(9j, is smooth; the difference is marginally 
significant with respect to the measurement accuracy. 



a level of accuracy and self-consistency compatible with the very demanding require- 
ments of the BepiColombo radioscience experiment. 

In particular, in this paper we have given the algorithm definitions for the compu- 
tation of the observables range and range rate, including the reference system effects 
and the Shapiro effect. We have shown which computation can be performed explic- 
itly and which ones need to be obtained from an iterative procedure. We have also 
shown how to push these computations, when implemented in a realistic computer 
with rounding-off, to the needed accuracy level, even without the cumbersome usage 
of quadruple precision. The list of "relativistic corrections", assuming we can distin- 
guish their effects separately, is long, and we have shown that many subtle effects are 
relevant to the required accuracy. However, in the end what is required is just to be 
fully consistent with a post-Newtonian formulation to some order, to be adjusted when 
necessary. Interestingly, the high accuracy of BepiColombo radio system may require 
implementation of the second post-Newtonian effects in range. 
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